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ABSTRACT 

A  hybrid  background  error  covariance  (BEC)  model  for  three-dimensional  variational  data  assimilation  of 
glider  data  into  the  Navy  Coastal  Ocean  Model  (NCOM)  is  introduced.  Similar  to  existing  atmospheric  hybrid 
BEC  models,  the  proposed  model  combines  low-rank  ensemble  covariances  Bm  with  the  heuristic  Gaussian- 
shaped  covariances  B0  to  estimate  forecast  error  statistics.  The  distinctive  features  of  the  proposed  BEC 
model  are  the  following:  (i)  formulation  in  terms  of  inverse  error  covariances,  (ii)  adaptive  determination  of 
the  rank  m  of  Bm  with  information  criterion  based  on  the  innovation  error  statistics,  (iii)  restriction  of  the 
heuristic  covariance  operator  B0  to  the  null  space  of  Bm,  and  (iv)  definition  of  the  BEC  magnitudes  through 
separate  analyses  of  the  innovation  error  statistics  in  the  state  space  and  the  null  space  of  B0. 

The  BEC  model  is  validated  by  assimilation  experiments  with  simulated  and  real  data  obtained  during 
a  glider  survey  of  the  Monterey  Bay  in  August  2003.  It  is  shown  that  the  proposed  hybrid  scheme  substantially 
improves  the  forecast  skill  of  the  heuristic  covariance  model. 


1.  Introduction 

In  recent  years,  development  of  hybrid  background 
error  covariance  (BEC)  models  has  been  an  area  of  active 
research  in  atmospheric  data  assimilation  (Hamill  and 
Snyder  2000;  Etherton  and  Bishop  2004;  Wang  et  al. 
2007).  It  has  been  shown  in  particular  that  hybrid  models 
tend  to  be  more  robust  than  conventional  ensemble- 
based  data  assimilation  schemes,  especially  when  the 
model  errors  are  larger  than  observational  ones  (Wang 
et  al.  2007,  2008,  2009).  This  feature  is  attractive  for  the 
regional  assimilation  problems  in  oceanography,  where 
information  on  the  background  state  is  often  scant  and 
incomplete. 

Sequential  data  assimilation  schemes  developed  so  far 
for  regional  oceanographic  studies  can  be  classified  in 
two  categories.  The  first  one  is  the  Kalman  filter  (KF)- 
type  algorithms  with  low-rank  BEC  matrices  Bm  derived 
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from  ensemble  statistics.  These  applications  encompass 
many  flavors  of  the  reduced-order  KF  techniques  (e.g., 
Evensen  2003;  Tippett  et  al.  2003;  Brasseur  and  Verron 
2006).  They  proved  to  be  especially  useful  for  monitor¬ 
ing  comparatively  large  domains  continuously  covered 
by  sea  surface  height/sea  surface  temperature  (SSH/ 
SST)  observations  at  the  surface  with  sporadic  vertical 
temperature/salinity  ( T/S )  soundings  by  Argo  drifters 
and  ships.  The  second  type  of  assimilation  algorithms  em¬ 
ploy  steady-state  covariances  B0  derived  from  long-term 
model  integrations  (Yin  et  al.  2011)  or  heuristic  Gaussian- 
shaped  covariance  operators  with  simple  dynamical  con¬ 
straints  (Weaver  and  Courtier  2001;  Pannekoucke  and 
Massart  2008).  The  latter  type  of  the  BEC  models  has 
recently  gained  considerable  attention  because  of  its  flex¬ 
ibility  and  convenience  in  introducing  prior  information 
into  the  covariance  model  in  cases  when  the  back¬ 
ground  model  solutions  are  biased  and/or  contain  large 
errors. 

A  typical  oceanographic  setting  of  such  kind  is  a 
near-coastal  survey  by  autonomous  gliders,  which  have 
recently  become  a  fast-developing  operational  technol¬ 
ogy  in  oceanography  (Rudnick  et  al.  2004).  Gliders  are 
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capable  of  making  remotely  controllable  surveys  of  lim¬ 
ited  areas  at  high  spatiotemporal  resolution.  Such  a  dense 
4D  coverage  is  usually  accompanied  by  a  relatively  poor 
knowledge  of  the  background  ocean  state:  near-coastal 
regions  are  often  affected  by  poorly  known  peculiari¬ 
ties  of  the  bottom  topography  and  the  associated  tidal/ 
inertial  motions  that  cannot  be  resolved  by  global  OGCMs. 
Considerable  model  error  covariances  also  persist  at  scales 
comparable  with  the  size  of  the  domain  due  to  incon¬ 
sistencies  in  the  boundary  conditions  and/or  local  atmo¬ 
spheric  forcing. 

Because  of  the  relative  novelty  of  glider  technology,  ex¬ 
amples  of  glider  assimilation  are  rare  in  literature  (Heaney 
et  al.  2007;  Shulman  et  al.  2009).  Recently,  Dobricic  et  al. 
(2010)  have  shown  that  three-dimensional  variational 
data  assimilation  (3DVar)  assimilation  of  glider  data 
significantly  improves  the  forecast  skill  of  a  regional 
model.  Most  importantly,  glider  data  were  able  to  cap¬ 
ture  basin-scale  BE  correlations,  which  improved  the 
model’s  forecast  skill  several  weeks  after  termination  of 
glider  observations.  Dobricic  et  al.  utilized  the  second 
category  3DVar  algorithm  based  on  stationary  Gaussian- 
shaped  BECs  in  the  horizontal  combined  with  EOF  de¬ 
composition  in  the  vertical  (Dobricic  and  Pinardi  2008) 
and  did  not  explicitly  include  adaptive  error  covariances 
inferred  from  model  statistics. 

In  this  study  we  propose  a  hybrid  3DVar  assimilation 
system  specifically  targeted  on  preserving  survey-scale 
correlations  that  could  be  resolved  by  gliders  in  coastal 
areas.  Similar  to  the  existing  atmospheric  hybrid  models, 
the  “flow  dependent”  part  of  the  covariance  Bm  is  defined 
as  a  low-rank  matrix  derived  from  ensemble  statistics. 
The  heuristic  part  of  the  covariance  B0  is  represented  by 
the  propagator  of  the  diffusion  equation  for  temperature 
and  salinity.  To  gain  extra  computational  efficiency,  the 
action  of  the  propagator  is  modeled  by  a  semi-implicit 
scheme  (Weaver  and  Ricchi  2004;  Yaremchuk  et  al.  2011, 
manuscript  submitted  to  Ocean  Modell.).  For  that  reason 
the  proposed  BEC  model  is  formulated  in  terms  of  the 
inverse  covariances  and  the  assimilation  problem  is 
solved  in  the  state  space  7 ZM . 

Another  distinctive  feature  of  the  BEC  model  is  an 
explicit  separation  of  the  covariance  components  in  7 ZM: 
the  action  of  B0  is  restricted  to  the  null  space  of  Bm.  This  is 
done  to  better  preserve  the  above-mentioned  regional- 
scale  error  correlations.  Since  low-rank  approximations 
of  large  covariance  matrices  tend  to  be  more  uncertain  at 
larger  distances  (Hamill  et  al.  2001),  we  paid  special  at¬ 
tention  to  the  determination  of  the  statistically  reliable 
number  of  modes  m  and  the  magnitudes  (scaling  co¬ 
efficients)  for  both  Bm  and  B0.  The  respective  algorithms 
are  based  upon  the  Bayesian  information  criterion  and 
analyses  of  the  innovation  statistics. 


The  rest  of  the  paper  is  organized  as  follows.  We  start 
with  the  description  of  the  hybrid  BEC  model  (section  2), 
then  briefly  review  the  Navy  Coastal  Ocean  Model 
(NCOM)  forecast  model  and  the  experimental  design 
for  the  Monterey  Bay  area  (section  3).  We  continue  with 
an  examination  of  the  forecast  skills  of  the  assimilation 
system  for  the  twin-data  experimental  setting  and  sub¬ 
sequent  real-data  experiment  (section  4).  Section  5 
concludes  the  paper. 

2.  A  hybrid  3DVar  assimilation  scheme 

a.  The  BEC  model 

The  analysis  increment  8x  of  the  sequential  data  as¬ 
similation  scheme  considered  here  is  obtained  by  mini¬ 
mizing  the  cost  function: 

J(8x )  =  ^[5xtB_15x+  (HSx  —  Sy)TR_1(HSx  —  8y)]  — >•  min, 

2  5x 

(i) 

where  B  is  the  BEC  matrix,  R  is  the  K  X  K  observation 
error  covariance  matrix,  T  denotes  transposition,  and  H 
is  the  linear  operator  projecting  model  state  x  E  7 ZM  on 
the  innovation  vector  8y  E  7 ZK,  whose  K  components  are 
the  model-data  misfits  of  the  background  solution. 

To  define  linear  operations  with  multivariate  vectors 
8x ,  we  introduce  a  diagonal  matrix  G  approximating  the 
background  error  variance.  Elements  of  G  depend  on  the 
physical  nature  of  the  fields  contributing  to  8x  and  spatial 
coordinates.  Farther  below  we  will  assume  that  all  the 
quantities  in  (1)  are  normalized  by  the  respective  error 
variances  and  introduce  new  variables: 

Sx*  =  G-1/25x,  8y*  =  R_1/2Sy. 

The  matrices  B,  H  are  appropriately  transformed  to 
keep  J  invariant: 

b;1  =  g'/2b  'gi/2,  h*  =  r1/2hg12. 

Dropping  the  asterisks  for  convenience  of  further  treat¬ 
ment,  the  cost  function  (1)  and  the  normal  equation 
8J/d8x  =  0  now  take  the  following  form: 

J(8x)  =  ^  [5xtB-1Sx  +  (H  8x  —  Sy)T(H5x  —  Sy)],  (2) 

[B_1  +  HtH]5x  =  HT8y.  (3) 

Farther  below  we  assume  R  to  be  known  and  focus  on 
the  structure  of  the  BEC  matrix  B. 

The  hybrid  covariance  models  developed  so  far  (Hamill 
and  Snyder  2000;  Etherton  and  Bishop  2004;  Wang  et  al. 
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2007)  utilize  linear  combinations  of  the  heuristic  or 
“static”  covariance  operator  B0  and  the  flow-dependent 
operator  Bm  derived  from  the  statistics  of  an  ensemble 
of  analyses/forecasts: 

B  =  «_lBm+r1B0.  (4) 


B  =  ^PAmPT  +  l[P±exp(— tD)PI]-\  (8) 

with  inversion  in  the  second  term  standing  for  the  gen¬ 
eralized  (Moore-Penrose)  inverse.  The  normal  equation 
(3)  takes  the  following  form: 


Here  a~1  and  /3_1  are  empirically  defined  positive  scalar 
parameters  often  constrained  by  the  requirement  a~x  + 
(3~1  =  1  (e.g.,  Wang  et  al.  2008).  We  adopt  the  tradi¬ 
tional  representation  of  Bm  in  the  form  of  a  matrix  de¬ 
fined  on  a  subspace  lZm  e  7 ZM  spanned  by  an  orthogonal 
basis  {e^}  derived  from  the  eigenvector  analysis  of  the 
ensemble  covariance: 


[aPAmxPT  +  /3P±  exp(— tD)PT  +  HTH]Sx  =  rfdy.  (9) 

Among  other  parameters,  the  BEC  model  (8)  depends 
on  the  inverse  magnitudes  a,  (3  of  its  components  and  the 
number  of  eigenvectors  m  spanning  the  lZm.  In  the  pro¬ 
posed  algorithm  a,  and  m  are  determined  from  the 
model  states  and  the  data. 


B 


m 


T 

PA  P  . 


m 


(5)  b.  Definition  of  m  and  a 


Here  P  is  a  rectangular  m  X  M  matrix  with  the  columns 
ek,  k  =  1,  . . .  ,  m  and  Am  is  the  diagonal  m  X  m  matrix 
whose  nonzero  elements  represent  the  variances  of  ek. 

In  the  absence  of  the  additional  prior  information,  B0 
is  often  represented  by  the  propagator  of  the  diffusion 
equation  (e.g.,  Weaver  and  Courtier  2001;  Pannekoucke 
and  Massart  2008): 

B0  =  exp(TD);  D  =  -VTv\,  (6) 

where  the  diffusion  tensor  v  depends  on  spatial  coordi¬ 
nates  to  simulate  inhomogeneity  and  anisotropy  of  the 
background  flow  and  r  is  the  scalar  parameter,  specifying 
in  3D  the  local  correlation  radii  pl  via  the  eigenvalues 
\lv  of  v:  pl  ~  \J2\1vt,  i  =  1,  2,  3.  The  parameter  t  can  be 
also  interpreted  as  “integration  time”  of  the  corresponding 
finite-difference  diffusion  equation. 

In  the  present  study,  we  adopt  the  diffusion  model  (6) 
and  define  the  inverse  of  the  BEC  operator  as 


Accurate  determination  of  the  first  term  in  the  BEC 
model  in  (8)  is  important  because  this  term  is  responsible 
for  capturing  error  correlations  on  scales  comparable 
with  the  size  of  the  domain.  In  oceanographic  applica¬ 
tions  these  errors  are  generated  by  poorly  known  open 
boundary  conditions  and  errors  in  atmospheric  forcing, 
which  tend  to  have  larger  scales  than  those  of  the  in¬ 
ternal  oceanic  variability.  In  addition,  Bm  may  contain 
valuable  information  on  the  dynamical  structure  of  the 
model  error  field  because  it  is  derived  from  the  prior 
statistics  of  the  forecast  errors. 

In  many  applications,  the  domain  surveyed  by  gliders 
is  rarely  well-observed  beforehand  and  the  first-guess  es¬ 
timate  of  the  background  state  may  be  far  from  reality.  So 
the  leading  eigenvectors  ek  of  the  first-guess  BEC  esti¬ 
mate  provide  poor  approximation  to  the  true  eigenvectors 
of  the  background  error  covariance.  To  assess  reliabil¬ 
ity  of  we  employ  the  Bayesian  information  criterion 
(Schwarz  1978)  and  define  the  optimal  number  of  “trus¬ 
ted”  eigenvectors  as  the  minimum  of 


B-^aPA-^  +  pP^Pl,  (7) 

where  P±  =  \M  —  PPT  is  the  projector  on  the  orthogonal 
supplement  of  lZm  and  \M  is  the  identity  operator  in  IZM . 
This  definition  statistically  separates  the  ensemble¬ 
generated  components  of  the  increment  PPT5x  from 
those  described  by  the  heuristic  BEC  model  B0.  Another 
reason  for  formulating  the  BEC  model  (7)  in  terms  of  the 
inverse  covariances  are  computational  advantages  of  the 
numerical  approximation  of  (6)  and  solving  the  normal 
equation  (3)  in  state  space  (Yaremchuk  et  al.  2011,  man¬ 
uscript  submitted  to  Ocean  Mo  dell). 

Since  B-1  has  a  two-cell  structure  in  an  orthogonal 
basis  containing  ek ,  the  respective  background  error 
covariance  matrix  can  be  readily  written  as 


C(m)  =  m  +  ^  In  o* 


min, 

m 


(10) 


where 


is  the  relative  residual  approximation  error  of  N  data 
samples  by  m  modes.  Here  rk  denotes  the  kth  obser¬ 
vation  location  at  the  nth  analysis  time,  the  coefficients 
/”  are  obtained  by  minimizing  the  numerator  in  (11)  at 
a  time  layer  n  for  a  given  number  of  modes  m  Kn,  and 
the  overbar  denotes  averaging  over  N  time  layers. 
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The  relationship  in  (10)  gives  an  asymptotic  (N  »  1) 
approximation  to  the  Bayesian  posterior  probability  for 
a  model  with  m  parameters  (linear  regression  on  m  ei¬ 
genfunctions)  given  N  observations  ( T/S  fields  sampled 
by  gliders  at  the  analysis  times)  under  the  assumption 
that  model-data  misfits  are  normally  distributed.  A 
similar,  but  less  restrictive  m  criteria  could  be  also  used 
(Akaike  1974;  Hannah  and  Quinn  1979). 

The  magnitude  of  Bm  is  determined  by  considering 
the  optimization  problem  (1)  in  the  ra-dimensional  sub¬ 
space  lZm  spanned  by  {e^}.  Because  the  Gaussian  part 
of  B  is  defined  in  the  orthogonal  supplement  of  7 Zm, 
an  approximate  formula  for  the  covariance  matrix  be¬ 
tween  the  projections  of  fix  on  ek  can  be  obtained  (see 
the  appendix): 


(SeSe1)  =  [«Am  +  Q] 


-1 


QA  Qt  +  Q 


[«A 


-1 


■or1. 

(12) 


Here  Q  =  PTHTHP  and  fie  is  the  m- dimensional  vector 
of  the  expansion  coefficients  such  that  fix  =  Pfie. 

Equation  (12)  can  be  used  to  compute  a  in  several 
ways.  The  matrices  Q  and  A  are  known  and  the  matrix 
on  the  left-hand  side  can  be  estimated  by  approximating 
fiy  by  the  linear  combinations  of  at  N  analysis  times  and 
computing  the  time-averaged  covariances  between  the 
vectors  fie  of  the  optimal  fit  coefficients.  Optimal  a  can  be 
then  computed  by  minimizing  a  norm  of  the  difference 
between  the  left-  and  right-hand  sides  of  (12).  Since  all  the 
matrices  in  (12)  are  positive  definite,  a  convenient  option 
is  to  set  the  difference  between  their  traces  to  zero.  In  the 
application  considered  below,  the  background  model  er¬ 
rors  are  much  larger  than  observational  errors  (|Q|  > 
a|Am-1|),  and  we  use  the  simplified  relationship 

(SeSeT)  ~  3  Am  (13) 

to  estimate  a:  its  value  is  found  by  minimizing  the  mean 
squared  difference  between  the  diagonal  elements  of 
(fiefieT)  and  Am/a. 

In  principle,  one  can  generalize  the  covariance  model 
in  1Zm  and  exactly  fit  the  observed  variances,  diag(fiefieT), 
by  adjusting  the  diagonal  elements  of  Am.  In  oceano¬ 
graphic  applications,  however,  there  is  no  reason  to  refine 
the  covariance  model  by  finetuning  the  eigenvalues  be¬ 
cause  even  the  leading  eigenvectors  of  Bm  are  known 
very  poorly.  Besides,  the  minimization  problem  is  non¬ 
linear  and  computationally  expensive.  We  therefore  choose 
a  simpler  model  (13)  with  a  single  scaling  factor  a. 


c.  Definition  of  f3 

Having  established  the  structure  of  the  dynamical  part 
of  the  covariance  model  (7)  we  can  now  determine  the 
magnitude  /3  of  the  Gaussian  part  by  equating  the  trace  of 
the  sample  forecast  error  covariance  Tr(fiyfiyT)  derived 
form  the  innovation  statistics  to  the  trace  of  HBHT  +  1^, 
a  technique  routinely  used  in  computation  of  the  inflation 
factor  in  the  Kalman  filtering  schemes  (e.g.,  Wang  et  al. 
2007).  Substituting  B  from  (8)  into  the  expression  HBHT, 
we  obtain 

(SyTSy)=Tr|3HPAmPTHT 

+  ^H[P±  exp(— tD)P^]_1Ht|  +  K, 

so  that 

_  Tr{H[P±  exp(— tD)P]J-1Ht j 
13  ~  <SyTSy>  -  K  -  Tr[HBmHT]/a  ' 


(14) 


(15) 


The  numerator  of  this  expression  can  be  computed  by 
the  Monte  Carlo  technique  (Bai  and  Golub  1997)  at  the 
expense  of  several  iterative  solutions  of  the  M  X  M 
system  of  equations  with  random  right-hand  sides. 

d.  Numerical  implementation 

In  the  present  study  we  used  a  simple  diagonal  model 
of  the  diffusion  tensor  v  =  diag  v\  assuming  that  local 
decorrelation  radii  pl  are  directly  proportional  to  the 
model  grid  steps  Ax',  i  =  1,  2,  3  spatially  varying  in  3D: 


Vvl  =  Ax*/V2.  (16) 

Since  gliders  directly  measure  only  the  temperature  and 
salinity  fields,  the  operator  B0  was  only  applied  to  the  T/S 
components  of  the  state  vector  under  the  prior  assump¬ 
tion  of  zero  correlations  between  them.  The  temperature/ 
salinity  background  error  correlations  were  taken  into 
the  account  by  the  Bm,  which  was  also  operating  in  the 
reduced  space  (i.e.,  eigenvectors  ek  were  only  estimated 
for  the  temperature  and  salinity  fields). 

The  BEC  operator  B0  in  (6)  was  approximated  by  an 
implicit  “time  integration  scheme”  (see,  e.g.,  Yaremchuk 
et  al.  2011,  manuscript  submitted  to  Ocean  Modell. ): 


exp(TD) 


M 


(17) 


where  r/n  is  the  length  of  the  implicit  “time  step”  and  n 
is  the  number  of  explicit  time  steps.  With  the  definition 
(16),  the  square  root  of  the  integration  time  t  has  the 
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Fig.  1.  2D  slices  of  the  rows  of  the  numerical  approximations 
of  exp(Dr).  The  patterns  were  obtained  by  applying  (a)  the  high- 
order  explicit  approximation  (I  +  Dt/lOO)100  and  (b)  the  implicit 
approximation  in  (16)  with  n  =  2  to  the  5-shaped  disturbances  of 
the  temperature  field  at  three  points  shown  by  white  circles.  The 
bathymetry  contours  are  in  m. 

meaning  of  the  proportionality  coefficient  between  the 
local  decorrelation  scales  pl  and  the  grid  steps  in  the  cor¬ 
responding  directions:  pl  =  yJ2 v^t  =  Axl^/r.  The  value 
of  t  —  20  was  determined  through  preliminary  experi¬ 
ments  described  in  section  3.  The  second  parameter  n  in 
(17)  was  chosen  as  a  compromise  between  the  numerical 
complexity  and  accuracy  in  approximation  of  exp(rD). 
For  n  =  2  the  approximation  error  is  close  to  15%  (Fig.  1), 


which  is  quite  reasonable  given  the  overall  uncertainty  in 
the  definition  of  the  heuristic  covariance  operator  B0.  At 
the  boundaries,  the  operator  D  was  specified  by  pre¬ 
scribing  zero  normal  derivatives. 

The  analysis  increment  Sx  was  obtained  by  solving  (9) 
with  a  generalized  minimum  residual  solver  (Saad  2003). 
In  correspondence  with  the  approximation  (17),  the  in¬ 
verse  of  exp(TD)  was  represented  by 


exp(rD) 


rD 

n 


(18) 


Depending  on  the  number  of  analyzed  observations  K , 
the  solution  of  (9)  required,  as  a  rule,  150-300  iterations, 
keeping  the  computational  cost  of  the  analysis  well  below 
the  cost  of  a  12-h  model  run  between  the  assimilations.  In 
the  reported  experiments,  the  state  space  dimension  M 
was  515  102,  which  is  the  number  of  observations  K 
varied  between  1500  and  3000,  and  m  never  exceeded  2. 
Within  these  ranges  of  K  and  m  the  CPU  time  required 
for  the  estimation  of  a  and  m  was  negligible  compared  to 
the  time  Tcpu  of  solving  the  normal  equation.  Conversely, 
estimation  of  p  required  several  (usually  4-7)  iterative 
inversions  of  P±  exp(-TD)Pl  at  the  expense  of  3-5tcpu, 
and  was  the  most  expensive  part  of  the  analysis. 

The  only  type  of  data  used  in  the  present  study  were 
temperature  and  salinity  profiles  from  gliders.  Therefore, 
balance  constraints  were  introduced  by  applying  the  lin¬ 
earized  equation  of  state  and  the  geostrophic-hydrostatic 
relationships  directly  to  the  temperature  and  salinity  in¬ 
crements  (e.g.,  Li  et  al.  2008)  obtained  from  minimization 
of  the  cost  function  (1). 


3.  Experiment  design 

The  BEC  model  was  verified  by  3DVar  assimilation  ex¬ 
periments  with  the  Navy  Coastal  Ocean  Model  (NCOM) 
configured  in  the  Monterey  Bay  (Fig.  2)  for  processing 
of  the  data  acquired  during  the  Autonomous  Ocean 
Sampling  Network  (AOSN  II)  experiment  (Ramp  et  al. 
2008).  The  experiment  was  conducted  in  the  summer  of 
2003  with  the  ultimate  goal  of  developing  an  adaptive 
sampling  technique  that  combines  numerical  forecasts 
with  the  data  flows  from  controllable  observation  plat¬ 
forms.  Observations  were  performed  by  several  types 
of  autonomous  underwater  vehicles  (AUVs)  including 
gliders,  high-frequency  radars,  two  moorings,  bottom- 
mounted  AD  CPs,  surface  drifters,  and  CTD  casts.  In  the 
present  study,  we  focus  the  analysis  on  the  temperature/ 
salinity  data  from  gliders  only:  space-time  coordinates  of 
the  gliders  are  used  to  define  observation  operators  H  (t) 
in  both  twin-  and  real-data  assimilation  experiments  with 
the  hybrid  3DVar  scheme. 
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Fig.  2.  Locations  of  glider  profiles  during  the  experiment  (solid 
dots)  and  model  grid  (smaller  dots).  The  bathymetry  contours  are 
in  m.  The  domain  used  for  estimation  of  the  distances  between 
the  model  states  in  twin-data  experiments  is  shown  by  the  solid 
black  line.  Circles  denote  locations  of  the  two  moorings  used  for 
validation  of  the  real-data  experiments. 


a.  Numerical  model,  observations,  and  validation 
technique 

To  simulate  oceanic  variability  during  the  experiment 
we  used  a  version  of  NCOM  forced  by  the  Coupled  Ocean- 
Atmosphere  Mesoscale  Prediction  System  (CO AMPS; 
Hodur  et  al.  2002)  winds  in  the  time  period  between 
1  August  ( t  =  0)  and  27  August  (t  =  27)  of  2003.  The 
model  was  configured  on  a  curvilinear  orthogonal  grid 
(Fig.  2)  with  horizontal  resolution  ranging  from  1  to 
4  km,  and  a  hybrid  crlz  vertical  coordinate  system  with  9 
cr  levels  in  the  upper  ocean  and  32  z  levels  below.  At  the 
open  boundaries,  the  model  was  one-way  coupled  to  the 
global  NCOM  model  (Shulman  et  al.  2009). 

Glider  observations  during  the  experiment  covered 
the  central  part  of  the  model  domain  (Fig.  2).  With  a  typ¬ 
ical  dive  cycle  of  about  1  h,  a  glider  would  travel  approx¬ 
imately  0.5  km  between  surfacings,  which  is  well  below 
the  grid  resolution.  For  that  reason  we  prescribed  obser¬ 
vational  operators  H  to  measure  instantaneous  vertical 
profiles  of  temperature  and  salinity  at  the  model  grid  point 
closest  to  the  average  of  the  surface  locations  of  a  glider 
before  and  after  a  dive.  In  the  assimilation  experiments  we 
used  a  12-h  analysis  cycle,  so  only  those  glider  profiles 
occurring  within  1-h  window  around  0000  and  1200  PST 
were  assimilated.  On  average,  the  model  domain  was 
covered  by  20-40  profiles  every  12  h. 


To  measure  distances  between  the  model  states,  a  di¬ 
agonal  metric  g  was  used.  The  diagonal  elements  of  g  ( gT , 
gs ,  gu,  gv,  and  gg)  were  depth  dependent  and  were  ob¬ 
tained  as  horizontally  averaged  time  variances  of  tem¬ 
perature  T,  salinity  S ,  horizontal  velocity  {u,  v},  and  SSH 
f ,  respectively,  at  a  grid  point  r: 

Sf(z)  =  <tf(r)-i(rj]2l\. 

In  the  above  equation,  £  stands  for  either  T,  S ,  u ,  v ,  or  £ 
and  angular  brackets  denote  the  horizontal  average  at 
level  z- 

Distances  rs  and  r8  between  the  model  states  were 
computed  in  both  observational  and  state  spaces: 

^(x1,x2)=((f1-f2)V)1/2; 

rj(xvx2)  =  {(€1-t2)2Rf)f.  (19) 

Here  the  angular  brackets  denote  averaging  over  the  3D 
model  domain  covered  by  gliders  (Fig.  2)  and  over  the 
glider  locations  rf ,  respectively. 

In  the  twin-data  experiments,  glider  “observations” 
of  temperature  yT  and  salinity  ys  were  extracted  from 
the  “true”  fields  T\  S*  (Fig.  3,  left  panel)  at  glider  lo¬ 
cations  rf  every  12  h  and  contaminated  by  white  noise  s 
with  zero  mean  and  0.1  rms  variation: 


yT  =  H  (T1  +  sgT );  ys  =  H(5'  +  egs). 

To  simulate  model  errors  and  assess  the  impact  of  as¬ 
similation  in  the  twin-data  experiments,  the  “first  guess” 
model  solution  xfg(r)  was  generated  by  integrating  the 
model  for  27  days  starting  from  the  initial  condition 
specified  by  x\t  =  8.5)  (Fig.  3,  right  panel). 

Using  xfg(0)  as  the  background  state  at  t  =  0,  a  series 
of  3DVar  assimilation  experiments  were  performed:  on 
every  12-h  assimilation  cycle,  model  forecasts  x^  were 
updated  with  the  analyses  increment  xa  =  xf  +  8x  and 
the  next  12-h  integration  was  started  from  xa.  The  skill 
of  assimilation  q(t)  was  assessed  in  both  observational 
and  state  spaces  by  calculating  the  normalized  distances 
between  the  12-h  model  forecasts  and  the  true  states: 


qf{t) 


rfiS^)  I, 

/fs(x',xfe)|0' 


(20) 


Experiments  with  real  data  were  conducted  and  the 
results  were  validated  in  a  similar  manner,  except  that  xf 
was  taken  as  the  first  guess  and  cfe  values  were  not  com¬ 
puted  because  the  true  state  was  unknown.  Instead  of  we 
assessed  the  forecast  skill  of  the  model  using  indepen¬ 
dent  temperature,  salinity,  and  velocity  observations  at 
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Fig.  3.  An  example  of  temperature  and  velocity  fields  for  the  (a) 
true  and  (b)  first-guess  solutions  (z  =  28  m,  t  =  1.5  days). 


two  moorings  shown  in  Fig.  2.  The  respective  distances 
rjp  and  skills  q ™  were  computed  similarly  to  (19)  and 
(20),  but  the  spatial  average  was  taken  over  two  moor¬ 
ings  and  the  differences  between  variables  in  (19)  were 
normalized  by  the  respective  rms  temporal  variations 
observed  at  the  moorings. 

b.  Parameters  of  the  hybrid  covariance  model 

In  contrast  to  atmospheric  applications,  regional  ocean¬ 
ographic  problems  have  more  difficulties  with  the  BEC 
estimation  from  ensembles.  The  reason  is  that  realistic 
ensembles  simulating  ocean  variability  on  regional  scales 
are  rarely  available.  In  the  present  study,  the  first -guess 
background  error  statistics  was  obtained  from  the  en¬ 
semble  of  the  differences  Sxf  between  the  first-guess 
states  x/g  at  times  enumerated  by  /  and  the  respective  12-h 


forecasts  (background  states)  x{  derived  from  the  assim¬ 
ilation  run  with  a  =  0.  The  latter  were  treated  as  a  rough 
approximation  to  the  ensemble  of  the  true  ocean  states. 
Since  the  number  of  ensemble  members  (1=1,...,  55) 
was  limited  by  the  duration  of  the  glider  survey,  the 
expected  number  of  statistically  sensible  eigenvectors  of 
the  ensemble  covariance  matrix  was  rather  small  and 
never  exceeded  2  (see  section  4). 

In  the  hybrid  assimilation  runs  (with  a  ^  0)  this  first- 
guess  ensemble  {8xf}  of  the  background  errors  was 
continuously  updated:  its  members  on  the  time  layers  / 
preceding  the  current  analysis  time  were  replaced  by  the 
members  derived  from  the  forecasts  already  made  with 
the  hybrid  scheme.  The  quality  of  updated  ensembles 
was  monitored  by  the  number  of  eigenvectors  accepted 
by  the  information  criterion  (10)  and  by  the  percentage 
of  8y  variance  explained  by  these  eigenvectors. 

To  increase  the  robustness  in  estimating  a  and  /3,  we 
utilized  the  method  of  Wang  et  al.  (2007)  and  performed 
additional  time  averaging  while  computing  the  sample 
variances  in  (13)  and  (15).  This  averaging  was  done  over 
the  ensemble  of  30  states  (15  days)  preceding  the  anal¬ 
ysis  time.  In  the  initial  15  days  of  the  assimilation  run, 
the  missing  background  states  were  taken  from  the  re¬ 
spective  forecasts  xfg(r)  generated  by  the  first-guess  so¬ 
lution.  Similar  averaging  over  N  =  30  samples  following 
the  analysis  time  was  done  when  estimating  crm  in  (11). 

There  are  two  parameters  in  the  definition  of  B0  that 
may  affect  the  performance  of  the  assimilation  scheme. 
One  is  the  “time  of  integration”  t  controlling  decorre¬ 
lation  length  scales,  and  the  other  is  the  order  n  of  ap¬ 
proximation  of  the  exponent  in  (17).  These  parameters 
were  tuned  by  numerical  experimentation. 

Comparison  of  the  model  solutions  (Fig.  3)  with  the 
grid  (Fig.  2)  gives  an  indication  that  horizontal  correla¬ 
tions  are  likely  to  decay  at  3-6  grid  steps.  We  checked 
this  hypothesis  by  twin  experiments  with  a  =  0  and 
computed  the  forecast  skill  of  the  assimilated  solutions 
with  various  values  of  y/r.  In  these  experiments,  n  was 
also  varied  in  the  range  between  1  and  4.  The  best  overall 
result  was  obtained  with  t  =  20  (pt  ~  4. 5 A xt)  and  n  =  2. 
Although  assimilation  quality  (with  r  =  15,  n  =  3)  was 
similar  and  in  some  periods  slightly  better,  the  compu¬ 
tational  cost  appeared  to  be  much  larger.  We  therefore 
selected  r  =  20,  and  n  =  2  as  basic  parameters  for  the 
assimilation  experiments. 

4.  Results 

a.  Twin-data  experiments 

Figure  4  compares  the  skill  of  3DVar  assimilation  runs 
performed  with  the  Gaussian  and  hybrid  BEC  models. 
During  the  first  8  days  of  assimilation,  the  hybrid  scheme 
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Fig.  4.  Normalized  distances  qss  and  qsv  between  the  true  solution  and  the  12-h  forecasts  of 
assimilated  solutions  with  a  =  0  (thin  lines)  and  a  /  0  (hybrid  model,  thick  lines).  The  plot 
below  shows  the  number  of  detected  eigenvectors  (bars,  right  axis)  and  the  portion  of  the 
model/data  misfit  variance  explained  by  those  modes  (thin  line  above  the  bars,  left  axis). 


was  unable  to  detect  any  statistically  reliable  modes. 
Between  days  8  and  11  the  first  mode  was  detected,  ac¬ 
counting  for  8%  of  the  forecast  error  on  day  8, 14%  on  day 
11  (Fig.  4),  and  17%  on  day  17.  On  day  12  the  second  mode 
was  detected,  accounting  for  4%  of  the  forecast  error  var¬ 
iance.  Contribution  of  the  second  mode  increased  to 
almost  10%  on  day  18.  Later,  the  modes  appear  to  lose 
their  predictive  skill  with  the  contributions  dropping  to 
12%  and  7%,  respectively,  on  day  25. 

The  12-h  forecast  errors  measured  in  terms  of  the 
normalized  distances  q \  from  the  true  state  are  found  to 
be  approximately  15%  smaller  than  for  the  assimilation 
run  with  a  =  0  (thin  lines  in  Fig.  4).  A  similar  level  of  error 
reduction  was  observed  by  Wang  et  al.  (2008)  in  twin-data 
experiments  with  a  hybrid  assimilation  into  the  Weather 
Research  and  Forecasting  (WRF)  model.  In  terms  of  the 
12-h  forecasts  of  temperature  and  salinity  in  glider  ob¬ 
servation  points  the  error  reduction  is  somewhat  smaller 
(11%  for  temperature  and  13%  for  salinity),  but  can  still 
be  considered  as  a  satisfactory  improvement  (Fig.  5). 


Assimilation  experiments  with  different  noise  in  ob¬ 
servations  have  shown  that  the  patterns  in  Figs.  4-5  are 
robust  up  to  the  noise  levels  of  0.5.  At  higher  noise  levels, 
the  approximation  (13)  becomes  less  accurate  and  it  is 
necessary  to  use  the  relationship  (12)  for  estimating  a. 
Larger  errors  in  estimating  a  result  in  the  loss  of  accuracy 
in  estimating  the  number  of  modes  m  and  the  magnitude  /3 
of  the  Gaussian  part  of  the  covariance.  We  therefore  as¬ 
sume  that  the  proposed  algorithm  is  valid  when  observa¬ 
tion  errors  are  considerably  smaller  than  the  background 
errors.  This  is  not  a  severe  restriction  for  regional  assim¬ 
ilation  problems  in  oceanography  where  the  first-guess/ 
background  model  solutions  are  rarely  preconditioned  by 
data  and  often  appear  to  be  rather  far  from  reality. 

b.  Real-data  experiments 

Figure  6  shows  a  typical  situation  we  encountered  in 
the  experiments  with  real  data  in  the  Monterey  Bay:  The 
first-guess  model  solution  does  not  have  much  in  common 
with  the  mooring  record  at  40  m  (left  panel).  Moreover, 


Fig.  5.  Normalized  distances  qgT  and  qss  between  the  true  solution  and  the  12-h  forecasts  of 
assimilated  solutions  with  a  =  0  (thin  lines)  and  a  ^  0  (hybrid  model,  thick  lines).  Distances  in 
observational  space  are  normalized  by  the  value  at  t  =  0.  The  thin  vertical  line  marks  the 
detection  time  of  the  first  mode  (see  Fig.  4). 
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Fig.  6.  (left)  Salinity  recorded  by  the  offshore  mooring  in  Fig.  2  (black  line)  and  the  corresponding  salinity  of  the  first-guess  NCOM 
solution  (gray  line),  (right)  Profiles  of  the  average  salinity  measured  by  gliders  (solid  bold  line),  moorings  (dashed  line),  and  extracted 
from  the  first-guess  model  solution  (solid  thin  line).  Shading  and  horizontal  bars  show  rms  variability. 


the  mean  profiles  in  the  right  panel  demonstrate  consid¬ 
erable  salinity  biases  above  30  m  and  in  the  depth  range 
between  50  and  200  m.  The  rms  variations  of  salinity 
measured  by  gliders  and  moorings  are  generally  consis¬ 
tent  with  each  other  in  magnitude  (cf.  horizontal  bars  and 
the  width  of  light  shading  around  the  thick  profile  in  the 
right  panel).  A  noticeable  bias  between  the  mean  sa¬ 
linity  measured  by  moorings  (solid  dots)  and  gliders 
(thick  line)  could  be  attributed  to  differences  in  aver¬ 
aging:  the  glider  profile  is  obtained  by  averaging  over  all 
the  glider  positions  (Fig.  2),  whereas  the  mooring  profile 
(solid  dots)  is  obtained  as  the  mean  of  only  two  moor¬ 
ings.  Similar  biases  between  the  first-guess  solution  and 
observations  were  obtained  for  the  temperature  field 
(not  shown). 

To  estimate  observation  errors,  we  compiled  the  glider 
T/S  records  at  times  when  gliders  passed  closer  than 
200  m  of  either  of  the  moorings  and  compared  these  data 
with  the  corresponding  observations  at  moorings.  In  to¬ 
tal,  168  of  such  “pass-by  events”  were  found.  Comparison 
of  these  observations  has  shown  that  the  rms  discrep¬ 
ancies  in  temperature  and  salinity  were  fairly  stable  with 
depth  and  varied  within  0.26-0.35  after  normalization  by 
the  rms  variances  am(z)  recorded  at  the  moorings  (hori¬ 
zontal  bars  in  the  right  panel  of  Fig.  6).  Based  on  these 
computations,  the  observation  error  variances  were  esti¬ 
mated  as  R1I2(z)  =  0.3crm(z)  and  assumed  not  to  vary  in 
the  horizontal. 

In  all  other  aspects  (the  first  guess,  the  background 
error  variance  G,  etc.),  assimilation  experiments  with 
real  data  were  configured  in  the  same  way  as  the  twin- 
data  experiments.  Because  the  “true  ocean  state”  in  the 
real-data  experiments  was  unavailable,  we  introduced 
an  additional  parameter  to  gauge  the  algorithm’s 
performance.  Similar  to  rf  was  computed  as  the 
normalized  distance  between  the  model  forecast  field  £ 


and  temperature,  salinity,  or  velocity  measured  at  the 
points  of  moored  observations: 

rz  =  ((gf-nW)m- 

Angular  brackets  denote  averaging  in  the  vertical  (11 
levels  of  T/S  observations  or  18  levels  of  AD  CP  data) 
and  over  two  moorings  shown  in  Fig.  2. 

Figure  7  demonstrates  the  differences  8q  between  the 
salinity  forecast  skills  q™ ,§  of  the  assimilation  run  with 
a  =  0  and  similar  forecast  skills  obtained  with  the  hybrid 
BEC  model.  Although  the  skill  improvement  does  not 
look  as  good  as  in  twin-data  experiments  (Fig.  4),  it 
appears  to  be  robust:  the  differences  in  skill  8q™  and  8q8s 
remain  positive  for  most  of  the  time  after  detection  of 
the  first  mode  on  day  3.  The  time  mean  values  for  8q 
8q™,  and  8q8T  were  found  to  be  1.3%,  2.8%,  and  2.0%, 
respectively. 

Compared  with  the  time-averaged  assimilation  skill  in 
twin-data  experiments  (e.g.,  qss  ~  0.4  for  salinity  in  Fig.  4), 
the  values  of  q^  in  real-data  experiments  were  much 
higher  (0.6-0.7  for  temperature/salinity  and  0.9  for  ve¬ 
locity).  This  difference  is  due  to  larger  observation  noise, 
its  more  complex  structure,  and  a  considerable  bias  (left 
panel  in  Fig.  6)  inconsistent  with  the  prior  statistical  as¬ 
sumptions.  The  Bayesian  algorithm  (10)  indicated  an 
occasional  presence  of  only  one  informative  mode:  de¬ 
tection  events  disappeared  on  day  14  and  reemerged  only 
at  the  end  of  the  assimilation  period  (Fig.  7).  Such  be¬ 
havior  could  be  attributed  to  the  poor  quality  of  the  first- 
guess  solution  and  insufficient  statistics  of  the  30-member 
ensemble  in  use.  Experiments  with  changing  the  ensem¬ 
ble  size  ne  (20-55  members)  have  shown  that  with  ne  = 
20,  the  number  of  detection  events  dropped  to  3,  whereas 
with  ne  =  50,  it  increased  from  30  (Fig.  7)  to  35  without 
any  substantial  improvement  of  the  forecast  skill.  Further 
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Fig.  7.  Improvement  of  the  12-h  salinity  forecasts  at  glider  observation  points  (gray  line)  and 
at  the  moorings  (black  line).  Positive  values  correspond  to  smaller  forecast  errors  for  the  hybrid 
scheme.  Vertical  bars  indicate  occasional  detection  of  only  one  mode  and  the  thick  black  line 
shows  the  percentage  of  the  error  variance  that  the  mode  explains. 


increase  of  the  ensemble  size  was  limited  by  the  duration 
of  the  assimilation  experiments. 

In  principle,  the  ensemble  could  be  expanded,  for  ex¬ 
ample,  by  the  breeding  technique,  but  the  problem  with 
the  poor  quality  of  the  first-guess  solution  (Fig.  6)  may 
still  persist,  because  the  bred  vectors  would  still  show 
unstable  modes  of  the  background  state  that  is  rather  far 
from  reality.  In  the  present  study,  we  used  a  simple  ap¬ 
proximation  to  the  error  fields  by  considering  an  en¬ 
semble  of  differences  between  a  free  model  run  and  an 
assimilation  run  with  the  Gaussian  covariance  model. 
This  ensemble  was  able  to  generate  just  a  few  members  in 
the  27-day  period.  One  may  hope,  however,  that  for 
longer  observation  periods  the  BEC  model  will  gain 
enough  skill  to  show  better  performance. 

Figure  8  shows  the  time  evolution  of  the  ratio  between 
the  weighting  parameters  a  and  p  in  the  twin-  and  real- 
data  experiments.  By  an  order  of  magnitude,  the  ratio  y  is 
consistent  with  the  results  of  Wang  et  al.  (2008)  who  set 
y  =  p/a  =  const  in  time  and  found  the  optimal  y  to  vary 
between  1  and  4  in  a  series  of  twin-data  experiments  with 
the  WRF  model.  In  our  case,  the  relative  weight  p  of  the 
Gaussian  term  in  the  cost  function  appeared  to  be  ap¬ 
proximately  2  times  smaller  in  the  twin-data  experiment 
(thin  curve  in  Fig.  7).  This  is  consistent  with  a  better  skill 
in  explaining  model-data  misfits  by  the  modes  retrieved 
in  the  twin-data  experiment  (cf.  Figs.  4  and  7).  Larger 
relative  values  of  a  on  days  10-13  (before  the  mode  re¬ 
jection)  can  be  explained  by  the  tendency  of  the  algo¬ 
rithm  to  keep  the  deteriorating  mode  “alive.” 

We  also  investigated  the  impact  of  the  algorithms  for 
definition  of  m  and  a  on  the  forecast  skill.  In  the  twin-data 
experiments  with  fixed  m  the  27-day-averaged  skill  was 
always  worse  than  that  in  Fig.  4  for  3  tested  values  of  y  = 
0.5,  1,  and  2.  When  m  was  computed  through  (10)  and 
y  was  kept  constant  at  0.95,  the  forecast  skill  was  virtually 
the  same  as  in  Fig.  4,  but  somewhat  below  using  other 


values  of  y.  Similar  results  were  obtained  with  real  data: 
keeping  m  =  const  degraded  the  forecast  skill,  often  be¬ 
low  the  one  obtained  with  Gaussian  BEC  model.  Several 
runs  with  an  adjustable  m  and  y  =  const  were  difficult  to 
interpret  as  the  skill  improvements  were  small,  highly 
variable,  and  did  not  show  any  deterministic  dependence 
on  the  value  of  y  e  [0.5,  2.5]. 

5.  Summary  and  discussion 

In  this  study  we  proposed  a  hybrid  BEC  model  spe¬ 
cifically  designed  for  3DVar  analysis  of  regional  circula¬ 
tions  supported  by  glider  surveys.  The  model  is  supplied 
by  an  algorithm  for  weighting  the  ensemble -generated 
error  covariance  Bm  against  the  heuristic  covariance  B0 
represented  by  the  propagator  of  the  diffusion  equation. 
Another  distinctive  feature  of  the  algorithm  is  the  de¬ 
tection  of  the  statistically  confident  eigenvectors  of  Bm  by 
means  of  the  Bayesian  information  criterion  (Schwarz 
1978).  The  method  is  based  on  the  assessment  of  the 
modes’  skill  in  approximation  of  the  forecast  error  fields 
accumulated  in  the  course  of  the  assimilation  run. 

The  proposed  BEC  model  is  formulated  in  terms  of  the 
inverse  covariances  with  the  restriction  of  B0  to  the  null 
space  of  Bm.  This  is  done  to  better  preserve  the  covari¬ 
ances  detected  by  the  information  criterion  and  captured 
by  Bm.  Formulation  of  the  minimization  problem  in  the 
state  space  allows  us  to  gain  extra  computational  effi¬ 
ciency  by  approximating  the  action  of  B0  via  a  semi- 
implicit  scheme. 

The  hybrid  BEC  model  was  validated  by  numerical 
experiments  with  simulated  and  real  data.  In  the  twin- 
data  setting,  the  hybrid  formulation  was  capable  of  im¬ 
proving  the  model’s  forecast  skill  by  15  %-20%,  which  is 
comparable  with  the  improvement  reported  by  Wang 
et  al.  (2008)  for  a  hybrid  scheme  with  the  atmospheric 
WRF  model.  Results  of  the  experiments  with  real  data 
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Fig.  8.  Time  variation  of  the  ratio  (3/a  in  the  twin-data  (thin  line)  and  real-data  experiments. 


showed  a  few  percent  improvement  with  sporadic  de¬ 
tection  of  only  one  mode.  We  attribute  this  to  a  poor 
quality  of  the  background  solution,  which  was  heavily  bi¬ 
ased  and  demonstrated  considerably  lower  time  variation 
in  the  temperature  and  salinity  fields  (left  panel  in  Fig.  6). 
Thus,  finding  a  better  background  solution  appears  to  be 
the  first  priority  in  upcoming  studies  of  the  algorithm. 

Other  developments  may  include  elaboration  of  the 
structure  of  the  diffusion  tensor  in  B0  and  improvement 
of  the  ensemble  generation  technique.  In  particular, 
diffusion  could  be  enhanced  along  the  f/H  contours  and / 
or  isopycnals  of  the  geostrophically  balanced  modes  if 
the  latter  are  detected.  The  ensemble  could  be  enriched 
by  the  vectors  bred  from  the  eigenvectors  of  Bm  or  just 
using  a  standard  breeding  technique,  if  the  background 
state  acquires  a  reasonable  forecast  skill  in  the  course  of 
assimilation.  Finally,  detected  eigenvectors  can  be  prop¬ 
agated  by  the  model  with  the  methods  used  in  Kalman 
filtering  schemes.  This  approach  will  naturally  combine 
the  advantages  of  statistical  and  dynamical  methods  and 
increase  versatility  of  the  hybrid  algorithm. 

One  of  the  drawbacks  of  the  proposed  model  is  the 
computational  cost  of  estimating  the  weight  (3  of  the  static 
covariance  [(15)].  Our  experience  shows,  however,  that 
the  numerator  in  (15)  weakly  depends  on  the  structure  of 
H  for  a  given  number  of  eigenvectors  and  can  be  effi¬ 
ciently  parameterized  by  a  linear  function  of  the  number 
of  observations.  In  fact,  the  predictive  skill  of  the  system 
did  not  change  when  such  a  linear  parameterization  was 
used.  A  similar  kind  of  parameterization  could  also  be 
employed  to  estimate  a  and  (3  when  the  background  er¬ 
rors  are  comparable  with  the  observation  errors. 

The  benefit  of  the  proposed  hybrid  model  may  also 
be  diminished  for  global  assimilation  problems  where 
some  sort  of  localization  is  needed  and  the  impact  of  the 
ensemble-generated  covariances  may  be  smaller  with 
higher  observation  density  (e.g.,  Whitaker  et  al.  2008). 
Nevertheless,  we  assume  that  the  proposed  approach  may 
have  a  prospect  for  further  development  for  regional  data 
assimilation  problems  with  poorly  known  background 
states. 
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APPENDIX 

Derivation  of  Equation  (12) 

Consider  the  optimization  problem  (1)  in  lZm  = 
span{e^}by  introducing  a  new  variable  8e  such  that  8x  = 
P8e: 

J  =  f  [5eTPTB-1PSe  +  (HPSe  -  Sy)T 

X  (HP5e  —  8y)]  — >  min  .  (Al) 

8eelZm 

The  normal  equation  (3)  is  now  reduced  to 

[PtB-1P  +  Q]8e  =  ESy,  (A2) 

where  the  operator  E  =  PTHT  projects  the  data  on  7 Zm 
and  Q  =  EET.  Substituting  the  adopted  inverse  co- 
variance  model  (7)  into  (A2)  and  taking  into  account  the 
identities  PTP  =  lm,  P±P  =  0,  the  normal  equation  is 
further  simplified  to 

[«Am  +  Q]'5e  =  ESy.  (A3) 

Introducing  the  notation  Y  =  (8y8yT)  for  8y  covariance, 
covariances  of  8e  could  be  written  as  follows: 

(SeSeT)  =  (aA"1  +  Q)“1EYET(aAm1  +  Q)_1.  (A4) 

On  the  other  hand,  in  accordance  with  the  observation 
model,  misfits  8y  between  the  background  state  and  the 
data  have  the  following  covariance: 

Y  =  HBHt  +  1^,  (A5) 

which,  after  projecting  on  the  eigenvectors  ek,  is 
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eyet  =  ehbht  +  eet 

=  ^AfflQT  +  ^EH[P1Bq1P1]_1HtEt  +  Q. 

(A6) 

Consider  now  the  middle  term  in  the  rhs  of  (A6): 

EH[P±Bo1Pir1HTET 

=  PtHtH[PB01PJ1HtHP  =  Q. 

For  “pointwise”  observations  (local  observational  op¬ 
erators)  the  matrix  HTH  is  equal  to  the  identity  matrix  \M 
with  diagonal  elements  masked  by  zeroes  in  the  points 
without  observations.  Therefore,  in  the  limit  of  a  per¬ 
fectly  observed  state  (HTH  =  lM),  this  term  vanishes.  For 
glider  observations,  which  densely  populate  the  domain 
during  the  survey,  one  may  assume  that  (HTH)  ~  \M  and 
neglect  this  term  in  the  time  average.  We  checked  the 
validity  of  this  assumption  by  estimating  the  ratio  |Q|/|Q| 
for  all  the  values  of  H  (t)  and  several  Ps  containing  the  first 
10  eigenmodes.  The  ratio  was  found  to  be  on  the  order  of 
10-2,  allowing  the  relationship  in  (A6)  to  be  reduced  to 

EYET=  lQAmQT  +  Q.  (A7) 


Substitution  of  (A7)  into  (A4)  yields  (12). 
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